* =====================================================
* Date: Dec 13, 2024
* Project: Health and anti-immigration attitudes
*
* Dataset: hlthmigr.dta
* 
*
* This file contains the code for the robustness tests.
* =====================================================

break

					* Prepare for analyses
					* ====================
**# Bookmark #1

*** Install packages
* ssc install coefplot
* ssc install cleanplots, replace
* ssc install psmatch2
* ssc install ivreg2
* ssc install ranktest

*** Define directory
clear all
set more off

* Insert your working directory here:
global PATH "..."

cd $PATH

use data_hlthmigr/hlthmigr.dta, clear
set scheme white_tableau

*** Define global variables for analyses			
** Demography
global Demogr	agea i.gndr2 i.mrtsts4 i.domicil3 i.brncntr2
** Socioeconomic factors 
global SES		eduyrs i.oesch5 hinctntb i.hincfel2 i.uempla i.hincsrc2a
** Political-cultural factors
global Cult		lrscale i.nwspol2 i.clsprty2 rlgdgr impsafe ipstrgv imptrad
** Clustering
global Clstr 	i.cntryc i.essround

					
					* Check for endogeneity
					* =====================


*** Propensity score matching

** Predict propensity score
probit hlth2a $Demogr $SES $Cult $Clstr [pweight=anweight]
est sto pscorew
predict pscorew

** Show common support, method #1
/* Kernel density graphs show decent evidence for common support */
#delimit	;
twoway 	(kdensity pscorew if hlth2a==0) (kdensity pscorew if hlth2a==1)
		, legend(order(1 "Untreated (Good health)" 2 "Treated (Bad health)") 
		pos(2) ring(0) region(fcolor(none)))
		ytitle(Kernel density) xtitle(Propensity score)
			;
#delimit cr

graph export output_hlthmigr/comsup1.emf, as(emf) replace

** Show common support, method #2
/* decent evidence for common support again. */
qui: psmatch2 hlth2a $Demogr $SES $Cult $Clstr, out(migrw100) common
psgraph

graph export output_hlthmigr/comsup2.emf, as(emf) replace

** Test coveratiate balance
* 1-1 matching w/o replacement
// Imbalance significantly reduced but some standardized difference remains

qui: psmatch2 hlth2a $Demogr $SES $Cult $Clstr, out(migrw100) common neighbor(1) noreplacement
pstest $Demogr $SES $Cult $Clstr, treated(hlth2a) both hist

* Multiple matches with replacement
// balance test standardized coefficient differences looking good, large differences have been reduced to differences below 10%; should ignore the significant differences based on t-tests (we have a pretty large sample with lots of not normally distributed covariates)

qui: psmatch2 hlth2a $Demogr $SES $Cult $Clstr, out(migrw100) common
pstest  $Demogr $SES $Cult $Clstr, treated(hlth2a) both hist

* Multiple matches with replacement and survey weights 
qui: psmatch2 hlth2a, pscore(pscorew) out(migrw100) common ate
pstest $Demogr $SES $Cult $Clstr, treated(hlth2a) both graph ndots(0)

** Estimate ATT 
* w/o survey weights 1-1 mathcing (1 neighbor) without replacement
psmatch2 hlth2a $Demogr $SES $Cult $Clstr, out(migrw100) common neighbor(1) noreplacement ate
matrix A = ( r(att) \ r(seatt) \ e(N) )

* w/o survey weights multiple matches with replacement 
psmatch2 hlth2a $Demogr $SES $Cult $Clstr, out(migrw100) common ate
matrix B = ( r(att) \ r(seatt) \ e(N) )

* with survey weights multiple matches with replacement 
psmatch2 hlth2a, pscore(pscorew) out(migrw100) common ate
matrix C = ( r(att) \ r(seatt) \ e(N) )

matrix define D = A,B,C
matrix colnames D = m1 m2 m3
matrix rownames D = "ATT" "SE" "No of observations"
estout matrix(D, fmt(%12.3fc)) using output_hlthmigr/psm.txt, replace

/* There is no straightforward way to automate exporting the results from psmatch2.
The easiest is to copy the treatment effect estimates into excel and format there. */



					* Robustness tests 
					* ================


*** Rerun main model in SEM with full information maximum likelihood (FIML)

/* The default FIML model using the mlmv option in Stata by default includes the 
means, variances, and covariances of observed exogenous variables among the 
parameters to be estimated. This approach is more appropriate acording to the 
Stata manual when handling large amout of missing data. However, this default 
FIML model is computationally too complex, the model does not converge. Therefore,
we will estimate the model without the means, variances, and covariances of the
exogenous variables. Using a smaller sample, the the two estimation methods yield 
identical results, so we conclude that it's appropriate to continue w/o estimating
the means, variances, and covariances of observed exogenous variables. */ 

#delimit	;
global 		Med1 agea gndr2 mrtsts4tab2-mrtsts4tab4 domicil3tab2 domicil3tab3 brncntr2
			eduyrs oesch5tab2-oesch5tab5 hinctntb hincfel2 uempla hincsrc2a
			lrscale nwspol2 clsprty2 rlgdgr impsafe ipstrgv imptrad
			cntryc_3-cntryc_27
			;
#delimit cr

cap drop imbgecor10 imuecltr10 imwbcntr10
gen imbgecor10=imbgecor*10 
gen imuecltr10=imuecltr*10
gen imwbcntr10=imwbcntr*10

#delimit	;
sem 		(Migrw -> imbgecor10 imuecltr10 imwbcntr10)
			(Migrw <- health $Med1 essround_2-essround_10)
			[pw=anweight]
			, method(mlmv) forcexconditional vce(cluster cntryround)
			;
#delimit cr
est sto robsem

#delimit ;
estout 	robsem
		using output_hlthmigr/robsem.txt 
		, replace cells((b(star fmt(%9.3fc)) se(fmt(%9.3fc)) z(fmt(%9.3fc)) p(fmt(%9.3fc)))) 
		stats(N, fmt(%9.0fc) labels("No. of observations" )) label 
		collabels("Coefficient" "S.E." "z" "p")
		eqlabels("Dep. var.: Worried about immigration")
		varlabels(
		_cons Constant  
		Migrw "Worried about immigration" 
		var(e.imbgecor) "var(e.Immigration is bad for the economy)"
		var(e.imuecltr) "var(e.Immigration is bad for culture)"
		var(e.imwbcntr) "var(e.Immigration makes country worse)"
		var(e.Migrw) "var(e.Worried about immigration)"
				) 
		legend mlabels(none) numbers
		refcat(agea "Demographic factors" eduyrs "Socioeconomic factors" 
		lrscale "Political-cultural factors" 
		oesch5tab2 "Social class (ref. cat.: Higher-grade service)"
		mrtsts4tab2 "Marital status (ref. cat.: Married or partnered)"
		domicil3tab2 "Domicile (ref. cat.: Big city or suburb)"
		var(e.trstpltr) "Variance components", nolabel)
		drop(cntryc* essround_*)
		;
#delimit cr

estat gof, stats(all) // get R-square (CD)

*** Robustness test with different dependent variable

** Replicate main OLS results with migrr100
* Not standardized, fully adjusted (b=1.397, p<0.001)
reg migrr100 health $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryround)

* Standardized coefficients

reg z_migrr100 z_health z_agea z_gndr2 z_mrtsts4tab2 z_mrtsts4tab3 z_mrtsts4tab4 ///
	z_domicil3tab2 z_domicil3tab3 z_brncntr2 ///
	z_eduyrs z_oesch5tab2-z_oesch5tab5 z_hinctntb z_hincfel2 z_uempla z_hincsrc2a ///
	z_lrscale z_nwspol2 z_clsprty2 z_rlgdgr z_impsafe z_ipstrgv z_imptrad ///
	$Clstr if in_pols [aw=anweight], cluster(cntryround)


** Test the mediation model using the three restict immigration items 

/* Four notes on these SEMSs:
	- We removed the ESS-round fixed effects which would make the models 
		close to impossible to estimate (testing such models after an hour still
		no sign of any convergence). The standardized coefficient in the unmediated 
		model is very close to the standardized OLS coefficient estimate, so this
		reassures us that removing the fixed effects does not introduce major bias.
		Structural equation models with so many dummies can be somewhat unstable,
		depending on your version of Stata and the exact parameters included, 
		you might have to experiment with the maximum set of country dummies 
		("cntryc_4-cntryc_26") to figure out what works. Unfortunately, these
		parameter estimates might show some slight variation depending on the 
		number of dummies included. The results reported in the article were 
		estimated using the full set of dummies, but this might not be viable in 
		all Stata versions or if the ESS files are updated.
	- There are two ways to calculate the latent variable.
		First, we can set the factor loading of one of the observed variables to 1. 
		Second, we can set the variance of the latent variable to 1. 
		This is the case when we use fully standardized model by specifying 
		the option standardized. This allows each of the (standardized) 
		factor loadings (i.e., correlations) to be calculated.
	- The default way of constructing the latent variable is listwise deletion.
		This is different from how alpha calculates a common index, which are
		used in the OLS models. Therefore, the sample sizes in the OLS and SEM
		models differ even though we are using the same set of variables. 
		This can be avoided by changing the estimation method to mlmv, which does
		not drop missing values. Combined with using in_pols as the estimation
		sample, we have exactly the same number of observations as in the OLS.
	- Regression table: there is no accepted way to present the regression table with
		SEM results, some fine-tuning of the estout table is needed.
*/

#delimit	;
sem			(Trstr -> trstpltr trstprlr trstlglr)
			(Migrr -> imsmetn imdfetn impcntr)
			(imblecor <- health $Med1)
			(ppltrstr <- health $Med1)
			(Trstr <- health $Med1)
			(Migrr <- health imblecor Trstr ppltrstr $Med1)
			if in_sem [pw=anweight]
			, standardized vce(cluster cntryround)
			;
#delimit cr
est sto sem4
estat eqgof

* Full regression table, appendix
#delimit ;
estout 	sem4
		using output_hlthmigr/sem4.txt
		, replace cells((b_std(star fmt(%9.3fc)) se(fmt(%9.3fc)) z(fmt(%9.3fc)) p(fmt(%9.3fc)))) 
		stats(N, fmt(%9.0fc) labels("No. of observations" )) label 
		collabels("Standardized coefficient" "S.E." "z" "p")
		eqlabels("Dep. var.: Immigrants drain public services" 
		"Dep. var.: Diminished interpersonal trust"
		"Dep. var.: Distrust in political institutions"
		"Dep. var.: Restrict immigration")
		varlabels(
		_cons Constant  
		Migrr "Restrict immigration" 
		Trstr "Distrust in political institutions"
		var(e.ppltrstr)  "var(e.Interpersonal trust)"
		var(e.trstpltr) "var(e.Trust in politicians)"
		var(e.trstprlr) "var(e.Trust in parliament)"
		var(e.trstlglr) "var(e.Trust in legal syststem)"
		var(e.imbgecor) "var(e.Immigration is bad for the economy)"
		var(e.imuecltr) "var(e.Immigration is bad for culture)"
		var(e.imwbcntr) "var(e.Immigration makes country worse)"
		var(e.imblecor) "var(e.Immigrants drain public services)"
		var(e.Trstr) "var(e.Distrust in political institutions)"
		var(e.Migrr) "var(e.Restrict immigration)"
				) 
		legend mlabels(none) numbers
		refcat(agea "Demographic factors" eduyrs "Socioeconomic factors" 
		lrscale "Political-cultural factors"
		oesch5tab2 "Social class (ref. cat.: Higher-grade service)"
		mrtsts4tab2 "Marital status (ref. cat.: Married or partnered)"
		domicil3tab2 "Domicile (ref. cat.: Big city or suburb)"
		var(e.trstpltr) "Variance components", nolabel)
		drop(cntryc_*)
		;
#delimit cr

* Standardized mediation path estimates (Table 3)
medsem, indep(health) med(imblecor) dep(Migrr) stand rit
medsem, indep(health) med(ppltrstr) dep(Migrr) stand rit
medsem, indep(health) med(Trstr) dep(Migrr) stand rit


*** Sensitivy to outliers
** Jacknife exclusion of each country

foreach i of num 1(1)30	{
	qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntryc!=`i' [aw=anweight], cluster(cntryround)
	est sto jacknife`i'
	}

#delimit ;
estout 	jacknife*
		using output_hlthmigr/jacknife.txt, replace 
		cells(b(star fmt("%9.3fc")) se(par fmt("%9.3fc")) p(par([ ]) fmt("%9.3fc"))) 
		stats(r2 N, fmt(%9.3fc %9.0fc)
		labels("R-Square" "No. of observations" )) label varlabels(_cons Constant) 
		legend nobaselevels noomitted numbers
		starlevels(* 0.05 ** 0.01 *** 0.001) 
		mgroups("Worried about immigration", pattern(1))
		indicate("Fully adjusted = *.essround *.cntryc agea *.gndr2 *.mrtsts4 
		*.domicil3 *.brncntr2 eduyrs *.oesch5 hinctntb *.hincfel2 *.uempla 
		*.hincsrc2a lrscale *.nwspol2 *.clsprty2 rlgdgr impsafe ipstrgv imptrad")
		refcat(1.hlth2a "Self-rated health (ref. cat.: Good)", nolabel)
		;
#delimit cr


** Exclude three countries with highest and lowest health
sort c_health_m
list c_health_m cntryc if atcntry
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="EE" & cntry!="LV" & cntry!="PT" [aw=anweight], cluster(cntryround)
est sto rob1
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="EE" & cntry!="LV" & cntry!="PT" [aw=anweight], cluster(cntryc)
est sto clust1

qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="IE" & cntry!="CY" & cntry!="GR" [aw=anweight], cluster(cntryround)
est sto rob2
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="IE" & cntry!="CY" & cntry!="GR" [aw=anweight], cluster(cntryc)
est sto clust2

** Exclude three countries with highest and lowest support for immigration
sort c_migr_m
list c_migr_m cntryc if atcntry
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="GR" & cntry!="CY" & cntry!="CZ" [aw=anweight], cluster(cntryround)
est sto rob3
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="GR" & cntry!="CY" & cntry!="CZ" [aw=anweight], cluster(cntryc)
est sto clust3

qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="LU" & cntry!="SE" & cntry!="IS" [aw=anweight], cluster(cntryround)
est sto rob4
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & cntry!="LU" & cntry!="SE" & cntry!="IS" [aw=anweight], cluster(cntryc)
est sto clust4

** Alternative measures of health
* Hampered by illness
qui reg migrw100 i.hlthhmp2 $Demogr $SES $Cult $Clstr [pweight=anweight] if in_pols, cluster(cntryround)
est sto rob5
qui reg migrw100 i.hlthhmp2 $Demogr $SES $Cult $Clstr [pweight=anweight] if in_pols, cluster(cntryc)
est sto clust5

* Unmet medical needs in the last year
qui reg migrw100 i.medtrunr $Demogr $SES $Cult i.cntryc [pweight=anweight], cluster(cntryc)
est sto rob6
est sto clust6

* Multiple health problems in the past year
qui reg migrw100 i.hltprm2 $Demogr $SES $Cult i.cntryc [pweight=anweight], cluster(cntryc)
est sto rob7
est sto clust7

* Unhealthy behaviors scale
qui reg migrw100 hlthbhv3 $Demogr $SES $Cult i.cntryc [pweight=anweight], cluster(cntryc)
est sto rob8
est sto clust8

reg migrw100 hlth2a happy $Demogr $SES $Cult i.cntryc [pweight=anweight], cluster(cntryround)
est sto rob9

** Robustness check with additional variables
* Family culture
qui reg migrw100 health i.medu i.fedu $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryround)
est sto rob9
qui reg migrw100 health i.medu i.fedu $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryc)
est sto clust9

* Childhood financial difficulties
qui reg migrw100 health i.emprf14 i.emprm14 i.fnsdfml3 $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryc)
est sto rob10
est sto clust10

* Labor market outsiderness (fixed-term contract)
qui reg migrw100 health i.wrkctra $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryround)
est sto rob11
qui reg migrw100 health i.wrkctra $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryc)
est sto clust11

* Economic insecurity
qui reg migrw100 health i.lknemny i.uemp5yr $Demogr $SES $Cult i.cntryc [aw=anweight], cluster(cntryround)
est sto rob12
qui reg migrw100 health i.lknemny i.uemp5yr $Demogr $SES $Cult i.cntryc [aw=anweight], cluster(cntryc)
est sto clust12

* Highest degree earned
qui reg migrw100 health i.edulevel $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryround)
est sto rob13
qui reg migrw100 health i.edulevel $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryc)
est sto clust13

* Restrict sample to native majority population
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & blgetmg2==0 & brncntr2==0 & migrfam==0 [aw=anweight], cluster(cntryround) 
est sto rob14
qui reg migrw100 health $Demogr $SES $Cult $Clstr if in_pols & blgetmg2==0 & brncntr2==0 & migrfam==0 [aw=anweight], cluster(cntryc)
est sto clust14 

* Another set of political attitudes variables
qui reg migrw100 health polintrr freehmsr $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryround)
est sto rob15
qui reg migrw100 health polintrr freehmsr $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryc)
est sto clust15

* Another set of conservative values variables
qui reg migrw100 health ipeqopt ipudrst ipcrtiv $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryround)
est sto rob16
qui reg migrw100 health ipeqopt ipudrst ipcrtiv $Demogr $SES $Cult $Clstr if in_pols [aw=anweight], cluster(cntryc)
est sto clust16

* Robustness check: country-round fixed effects
qui reg migrw100 health $Demogr $SES $Cult $Clstr i.cntryround if in_pols [aw=anweight], cluster(cntryround)
est sto rob17
qui reg migrw100 health $Demogr $SES $Cult $Clstr i.cntryround if in_pols [aw=anweight], cluster(cntryc)
est sto clust17

* Robustness check: without survey weights
qui reg migrw100 health $Demogr $SES $Cult $Clstr i.cntryround if in_pols, cluster(cntryround)
est sto rob18
qui reg migrw100 health $Demogr $SES $Cult $Clstr i.cntryround if in_pols, cluster(cntryc)
est sto clust18

cap est drop robsem

#delimit ;
estout 	rob*
		using output_hlthmigr/rob.txt, replace 
		cells(b(fmt("%9.3fc")) se(par fmt("%9.3fc")) p(par([ ]) fmt("%9.3fc"))) 
		stats(r2 N, fmt(%9.3fc %9.0fc)
		labels("R-Square" "No. of observations" )) label varlabels(_cons Constant) 
		legend nobaselevels noomitted numbers
		mgroups("Worried about immigration", pattern(1))
		keep(health *.hlthhmp2 *.medtrunr *.hltprm2 hlthbhv3) 
		refcat(1.hlthhmp2 "Hampered by illness (ref. cat.: No)"
		1.medtrunr "Unmet medical needs last year (ref. cat.: No)"
		1.hltprm2 "Multiple health problems last year (ref.cat.: No)", nolabel)
		;
#delimit cr


#delimit ;
estout 	clust*
		using output_hlthmigr/clust.txt, replace 
		cells(b(fmt("%9.3fc")) se(par fmt("%9.3fc")) p(par([ ]) fmt("%9.3fc"))) 
		stats(N_clust N, fmt(%9.0fc %9.0fc)
		labels("No. of countries" "No. of observations" )) label varlabels(_cons Constant) 
		legend nobaselevels noomitted numbers
		mgroups("Worried about immigration", pattern(1))
		keep(_cons)
		;
#delimit cr

*** Checking the robustness of the interaction effect
** Use national-level perceived share of immigrants (noimbro)

* Split noimbro into top and bottom quartiles
xtile noimbro_quartile = noimbro, nq(4)
gen noimbro_topbottom = .
replace noimbro_topbottom = 1 if noimbro_quartile == 1 // Bottom 25%
replace noimbro_topbottom = 2 if noimbro_quartile == 4 // Top 25%
cap lab drop noimbro_tb_lbl 
label define noimbro_tb_lbl 1 "Bottom 25%" 2 "Top 25%"
label values noimbro_topbottom noimbro_tb_lbl


* Plot interaction top and bottom 25%
reg migrw100 c.health##noimbro_topbottom $Demogr $SES $Cult $Clstr, cluster(cntryround)
est sto area
margins, at(health=(1(1)5) noimbro_topbottom=(1 2))

#delimit 		;
marginsplot, 	plot1opts(lcolor(navy%90) mcolor(navy%90)) ci1opts(lcolor(navy%90))
				plot2opts(lcolor(cranberry%90) mcolor(cranberry%90)) ci2opts(lcolor(cranberry%90))
				title("") ytitle("Worried about immigration") xlabel(,angle(45) labsize(*0.9)) 
				xtitle("Self-rated health") 
				xsc(titlegap(*-2))
				legend(cols(1) rows(5) pos(10) ring(0) size(small)
				bmargin(zero) rowgap(0) colgap(0) region(fcolor(none)) 
				title("Selected values of the moderator" 
				"(Perceived national prevalence of immigrants):", 
				j(left) size(medsmall)))
				;
#delimit cr
graph export output_hlthmigr/interactionrob.emf, as(emf) replace


					* Test the validity of composite variables
					* ========================================

					
** Test if main mediator and proposed dependent variable are distinct enough
factor imuecltr imbgecor imwbcntr, ml

esttab using output_hlthmigr/factor1.rtf, cells("L[1](t fmt(%9.3fc) label(Factor 1)) Psi(fmt(%9.3fc) label(Uniqueness))") nogap noobs nonumber nomtitle label rtf replace

factor imuecltr imbgecor imwbcntr imblecor, ml

esttab using output_hlthmigr/factor2.rtf, cells("L[1](t fmt(%9.3fc) label(Factor 1)) Psi(fmt(%9.3fc) label(Uniqueness))") nogap noobs nonumber nomtitle label rtf replace

** sem produces same results as factor with ml estimation method
/* note: our main models including the structural components differ slightly
because there we cluster standard errors and use a sample restricted to variables
with full information (this might change if I change to fiml).*/

* Set shared sample
qui sem (Migrw -> imbgecor imuecltr imwbcntr imblecor)

cap drop in_migrw
gen in_migrw = e(sample)

* Factor model 1
sem (Migrw -> imbgecor imuecltr imwbcntr) if in_migrw, standardized
est sto migrw1

#delimit ;
esttab 	migrw1
		using output_hlthmigr/migrw1.rtf
		, replace rtf cells((b_std(fmt(%9.3fc)))) 
		stats(N, fmt(%9.0fc) labels("No. of observations" )) label 
		collabels("Factor loading")
		varlabels(Migrw "Worried about immigration") 
		legend mlabels(none) numbers drop(_cons var(e.imbgecor) var(e.imuecltr) 
		var(e.imwbcntr) var(Migrw))
		;
#delimit cr

estat eqgof
estat gof, stats(all)

*Factor model 2
sem (Migrw -> imbgecor imuecltr imwbcntr imblecor) if in_migrw, standardized
est sto migrw2

#delimit ;
esttab 	migrw2
		using output_hlthmigr/migrw2.rtf
		, replace rtf cells((b_std(fmt(%9.3fc)))) 
		stats(N, fmt(%9.0fc) labels("No. of observations" )) label 
		collabels("Factor loading")
		varlabels(Migrw "Worried about immigration") 
		legend mlabels(none) numbers drop(_cons var(e.imbgecor) var(e.imuecltr) 
		var(e.imwbcntr) var(e.imblecor) var(Migrw))
		;
#delimit cr

estat eqgof
estat gof, stats(all)
